#include <stdio.h>
#include <gmp.h>
#include <string.h>
#include <sys/time.h>
#include <unistd.h>
#include <pthread.h>

//DEFINES

//BOOLEANAS PARA O TESTE DE CONDICOES
#define FALSE  0
#define TRUE 1



typedef struct{
    mpf_t Y[2];
    mpf_t A[2];
    mpf_t PI[2];
    int i;
    int conv;
}termo_borwein;

//PRECISAO = 0,000000000 .... 0000000001 (GLOBAL)
mpf_t prec;



/* Verifica a convergência de dois valores utilizando a variável prec */
//SE PI_2 - PI_1 < PRECISAO, RETORNA 1 (CONVERGIU) SENAO RETORNA 0 (NAO ESTA NA PRECISAO
int verifica_convergergencia(mpf_t PI_1, mpf_t PI_2)
{
    mpf_t aux;
    int resposta;

    mpf_init(aux);

    // aux = |PI_2 - PI_1|
    mpf_sub(aux, PI_2, PI_1); mpf_abs(aux, aux);

    //COMPARA A PRECISAO
    resposta = (mpf_cmp(aux, prec) <= 0);

    mpf_clear(aux);

    return resposta;
}




//FUNCAO QUE RODARA EM UMA THREAD, PARA O CALCULO DO FATOR Y DA EQUACAO DE BORWEIN
void *y_paralelo(void *termo)
{
    termo_borwein *bo = (termo_borwein*) termo;

    //PARA SABER SE O QUAL O TERMO MAIS ANTIGO, USA-SE O MODULO DO VALOR DA ITERACAO ATUAL
    int i0 = (bo->i) % 2;
    int i1 = (bo->i + 1) % 2;

    mpf_t aux1, aux2;
    mpf_inits(aux1, aux2, NULL);

    //aux1 = (1-  y0^4)^(1/4)
    mpf_pow_ui(aux1, bo->Y[i0], 4);
    mpf_ui_sub(aux1, 1, aux1);
    mpf_sqrt(aux1, aux1);
    mpf_sqrt(aux1, aux1);

    //aux2 = aux1 + 1
    mpf_add_ui(aux2, aux1, 1);

    //y1 = (1 - aux1)/aux2
    mpf_ui_sub(bo->Y[i1], 1, aux1);
    mpf_div(bo->Y[i1], bo->Y[i1], aux2);


    //LIMPA VARIAVEIS GMP
    mpf_clears(aux1, aux2, NULL);


    return NULL;
}

//FUNCAO QUE RODARA EM UMA THREAD, PARA O CALCULO DO FATOR A DA EQUACAO DE BORWEIN
void *a_paralelo(void *termo)
{
    termo_borwein *bo = (termo_borwein *) termo;

    //SE FOR A PRIMEIRA ITERACAO, DEVE-SE ESPERAR O PRIMEIRO TERMO Y
    if(bo->i < 1){
        return NULL;
    }

    //PARA SABER A ORDEM CICLICA DOS FATORES
    int i1 = (bo->i-1 + 1) % 2;
    int j0 = (bo->i-1) % 2;
    int j1 = (bo->i-1 + 1) % 2;


    //CALCULO DE N
    int n = 2*(bo->i-1)+3;


    //VARIAVEIS AUXILIARES
    mpf_t aux3, aux4;
    mpf_inits(aux3, aux4, NULL);

    //aux1 = y1*(1+y1 + y1^2)
    mpf_pow_ui(aux3, bo->Y[i1], 2);
    mpf_add(aux3, aux3, bo->Y[i1]);
    mpf_add_ui(aux3, aux3, 1);
    mpf_mul(aux3, aux3, bo->Y[i1]);

    //aux2 = aux1*2^n
    mpf_set_ui(aux4, 2);
    mpf_pow_ui(aux4, aux4, n);
    mpf_mul(aux4, aux4, aux3);

    //a1 = a0*(1+y1)^4 - aux2
    mpf_add_ui(bo->A[j1], bo->Y[i1], 1);
    mpf_pow_ui(bo->A[j1], bo->A[j1], 4);
    mpf_mul(bo->A[j1], bo->A[j1], bo->A[j0]);
    mpf_sub(bo->A[j1], bo->A[j1], aux4);




    //LIMPEZA DAS VARIAVEIS
    mpf_clears(aux3, aux4, NULL);



    return NULL;
}

//FUNCAO QUE RODARA EM UMA THREAD, PARA O CALCULO DO PI
void *pi_paralelo(void *termo)
{
    termo_borwein *bo = (termo_borwein *) termo;

   //SE FOR A PRIMEIRA OU SEGUNDA ITERACAO, DEVE-SE ESPERAR PARA COMPLETAR O PIPELINE
    if(bo->i < 2){
        return NULL;
    }

    int j0 = (bo->i-2) % 2;
    int j1 = (bo->i-2 + 1) % 2;

    //pi = 1/a
    mpf_ui_div(bo->PI[j1], 1, bo->A[j1]);




    //CASO HAJA CONVERGENCIA, CONV EH SETADA COMO TRUE
    if(verifica_convergergencia(bo->PI[j0], bo->PI[j1])) bo->conv = TRUE;


    return NULL;
}



//METODO DE BORWEIN CONCORRENTE
void borwein_concorrente(mpf_t pi, long *it)
{
    //THREADS
    pthread_t Thread[3];

    //TIPO DE DADOS COM VARIAVEIS DE BORWEIN
    termo_borwein bo;

    //AUXILIARES
    int i, j;



    //INICIALIZACAO DE VARIAVEIS COM OS TERMOS DE BORWEIN
     // A0 = 6 - 4 * sqrt(2)
    mpf_init_set_ui(bo.A[0], 2);
    mpf_sqrt(bo.A[0], bo.A[0]);
    mpf_mul_ui(bo.A[0], bo.A[0], 4);
    mpf_ui_sub(bo.A[0], 6, bo.A[0]);

    // Y0 = sqrt(2) - 1
    mpf_init_set_ui(bo.Y[0], 2);
    mpf_sqrt(bo.Y[0], bo.Y[0]);
    mpf_sub_ui(bo.Y[0], bo.Y[0], 1);

    //INICIA AS VARIAVEIS GMP
    mpf_inits(bo.Y[1], bo.A[1], bo.PI[0], bo.PI[1], NULL);

    //SE A O PI CONVERGIU PARA O VALOR COM A PRECISAO DEFINIDA
    bo.conv = FALSE;


    //ALGORITMO ITERATIVO (ENQUANTO NAO ATINGIR A CONVERSAO DEFINIDA)
    for(i = 0; !bo.conv; i++)
    {
        //NUMERO DA ITERACAO
        bo.i = i;

        //CRIA UMA THREAD PARA CADA PARAMETRO A, Y, PI
        pthread_create(&Thread[0], NULL, y_paralelo, &bo);
        pthread_create(&Thread[1], NULL, a_paralelo, &bo);
        pthread_create(&Thread[2], NULL, pi_paralelo, &bo);

        //SINCRONISMO DAS THREADS (ESPERA TODAS ACABAR)
        for(j = 0; j < 3; j++)
            pthread_join(Thread[j], NULL);
    }


     //PASSA O NUMERO DE ITERACOES POR REFERENCIA
    *it = i;

    //SETA O NOVO VALOR DE PI
    mpf_set(pi, bo.PI[(i + 1)%2]);

    //LIMPA VARIAVEIS GMP
    mpf_clears(bo.Y[0], bo.Y[1], bo.A[0], bo.A[1], bo.PI[0], bo.PI[1], NULL);
}




int main(int argc, char *argv[])
{
    mpf_t pi;
    long it;
    int time;

    //VARIAVEIS PARA GUARDAR O TEMPO
    struct timeval inicio, fim;

    //SETAR A PRECISAO DO NUMERO DE CASAS DO TIPO DE DADOS MPG
    //EM BITS = (NUM_DIGITS * log2(10) ) = 33219281
    //ISSO EH NECESSARIO, POIS A PRECISAO DA FUNCAO MPG EH DADA EM VALOR DE BITS

    mpf_set_default_prec(33219281);
    // prec = 1/(10^PREC)
    mpf_init_set_ui(prec, 10);
    mpf_pow_ui(prec, prec, 33219281); mpf_ui_div(prec, 1, prec);
    //INICIA VARIAVEL PI
    mpf_init(pi);


    printf("RODANDO METODO DE BORWEIN CONCORRENTE...\n");

    //PEGA O TEMPO INICIAL DE EXECUCAO
    gettimeofday(&inicio, NULL);


    //METODO PRINCIPAL
    borwein_concorrente(pi, &it);


    //PEGA O TEMPO FINAL DE EXECUCAO
    gettimeofday(&fim, NULL);

    //DIFERENCA ENTRE FINAL E INICIAL
    time = fim.tv_sec - inicio.tv_sec;

    //SALVA O PI EM UM ARQUIVO TXT
    FILE *fp = fopen("borwein_concorrente.txt", "w");
    // IMPREIME COM PRECISAO  DE 10000000 DIGITOS*/
    mpf_out_str(fp, 10, 10000000, pi);
    fclose(fp);



    printf("COMPLETO\nITERACOES: %ld\n", it);
    printf("TEMPO DE EXECUCAO: %d segundos\n", time);

    //LIMPEZA DE VARIAVEIS GMP
    mpf_clear(prec);
    mpf_clear(pi);

    return 0;
}
